/****************************************************************************************************
* Title: Workhorses of Opportunity, Howard and Weinstein
*
* Inputs:
*   - justnormasylum.dta
*   - county_outcomes_dta.dta
*
* Outputs:
*   - xxx_with_average_by_sex.pdf (sex-disaggregated graphs of normal school effects)
*
* Description:
*   Estimates the effect of being from a county with a normal school on various outcomes
*   across income percentiles, separately by sex. Plots these effects alongside baseline values.
****************************************************************************************************/

* Load and prepare normal school indicator
clear 
use "justnormasylum"
rename hasnormalschool hasnormal
keep hasnormal cty_fips 
tempfile normal
save `normal'

* Set up environment
clear
clear matrix
clear mata
set maxvar 15000

* Load Chetty data and merge with normal school flag
clear
use county_outcomes_dta
gen cty_fips = state * 1000 + county
rename state statefips
rename county countychetty
merge 1:m cty_fips using `normal', keep(3)

* Initialize outcome storage
gen outcome = .
gen outcome0 = .
tempfile results

* Estimate treatment effects by sex and income percentile
foreach xxxxx in has_dad has_mom hs coll comcoll somecoll married jail staycz stayhome working kfr kir {
	foreach percentile in 1 25 50 75 100 {
		foreach sex in male female {
			
			replace outcome = `xxxxx'_pooled_`sex'_p`percentile'
			replace outcome = outcome * 100 if inlist("`xxxxx'", "kfr", "kir")
			
			reghdfe outcome 1.hasnormal, absorb(statefips) cluster(statefips) nocon
			
			summ outcome if !hasnormal
			local basemean = r(mean)
			
			preserve
				parmest, fast
				gen percentile = `percentile'
				gen sex = "`sex'"
				gen outcome = "`xxxxx'"
				gen baselinemean = `basemean'
				cap append using `results'
				save `results', replace
			restore
		}
	}
}

* Prepare data for plotting
clear
use `results'
gen min90 = -invttail(dof, .05) * stderr + estimate
gen max90 =  invttail(dof, .05) * stderr + estimate

* Generate plots by sex for each outcome
foreach xxxxx in hs coll somecoll married jail staycz stayhome working kfr kir {
	preserve
		keep if outcome == "`xxxxx'"
		replace percentile = percentile + 1 if sex == "female"

		twoway rspike min95 max95 percentile if sex == "female" || ///
			rcap min90 max90 percentile if sex == "female", color(dkgreen) || ///
			scatter estimate percentile if sex == "female", mcolor(dkgreen) yline(0, lpattern(dash) lcolor(gs8)) ///
			ytitle("Effect of Normal County (Spikes)") xtitle("Parents' Income Percentile") xlabel(0 25 50 75 100) || ///
			rspike min95 max95 percentile if sex == "male", color(maroon) || ///
			rcap min90 max90 percentile if sex == "male", color(maroon) || ///
			scatter estimate percentile if sex == "male", mcolor(maroon) || ///
			line baselinemean percentile if sex == "female", lcolor(dkgreen*.5) lpattern(dash) yaxis(2) ///
			ytitle("Average Value in Asylum County (Dashed)", axis(2)) name(`xxxxx'_average, replace) || ///
			line baselinemean percentile if sex == "male", yaxis(2) lcolor(maroon*.5) lpattern(dash) ///
			legend(order(3 6) label(3 "Female") label(6 "Male"))

		graph export "`xxxxx'_with_average_by_sex.pdf", as(pdf) replace
	restore
}
